Energy Analytics (Online) Risk and Hedging Assignment: 2020

Longin Le - CID (01626009)

In [100]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import norm
import random
from scipy.stats import multivariate_normal
In [101]:
data = pd.read_excel("Data for risk and hedging assignment.xlsx", index_col=0, names=['Day','gas','electricity'])

# Changing units from pences to pounds 
data.gas = data.gas/100
In [102]:
data.head()
Out[102]:
gas electricity
Day
1 0.430696 53.71
2 0.430608 47.69
3 0.437765 54.50
4 0.434990 58.77
5 0.437415 50.00

Part 1

In this part I will estimate the average daily profit over the 30-day period, based on the market prices for gas and electricity from the given dataset for November 2019.

Assumptions:

  • Drango operates from 8AM to 8PM - 12 hours per day
  • Drango consumes 48'222.65 gas therms per hour, thus 12 hours * 48'222.65 therms/hour = 578'671.80 gas therms per day at full capacity
  • A therm of gas will produce 14.516 kWh of energy
  • Drango produces maximum capacity of energy of 12 hours 48'222.65 gas therms/hour 14.516 kWh/ gas therm = 8'399.99 MWh a day
  • Drango gets fixed amount of 250'000 therms of gas from the fixed price contract at 0.5 GBP per day
  • Gas therms is purchased at market prices - 578'671.89 - 250'000 = 328'671.80 gas therms per day
  • Drango daily fixed operation costs are 140'000 GBP

Calculations:

  • Average daily profit = profit for a day / number of days
  • Profit for a day = Sales - Total Cost
  • Sales = 8'400 MWh * market price for a given day
  • Total Cost = Gas consumption with contract price + Gas consumption with market price + Operation daily
  • Gas consumption with contract price = 250'000 therms * 0.5 GBP = 125'000 GBP / day
  • Gas consumption with market price = 328'671.80 therms * market price for a given day
  • Operation daily cost = 140'000 GBP
In [103]:
# Copy of the dataset 
df_part1 = data.copy()
In [104]:
# Gas consumption with contract price
df_part1['Gas cost (contract)'] = 125000

# Gas consumption with market price
df_part1['Gas cost (market)'] = (578671.80-250000)*df_part1.gas

# Operation daily cost
df_part1['Operation cost'] = 140000

# Total Cost
df_part1['Total Cost'] = df_part1['Gas cost (contract)'] + df_part1['Gas cost (market)'] + df_part1['Operation cost']

# Sales
df_part1['Sales'] = (8399999.85/1000)*df_part1.electricity

# Profit
df_part1['Profit'] = df_part1['Sales']-df_part1['Total Cost']

df_part1.round(2).style.format("{:,}")
Out[104]:
gas electricity Gas cost (contract) Gas cost (market) Operation cost Total Cost Sales Profit
Day
1 0.43 53.71 125,000 141,557.59 140,000 406,557.59 451,163.99 44,606.4
2 0.43 47.69 125,000 141,528.79 140,000 406,528.79 400,595.99 -5,932.8
3 0.44 54.5 125,000 143,881.04 140,000 408,881.04 457,799.99 48,918.95
4 0.43 58.77 125,000 142,968.95 140,000 407,968.95 493,667.99 85,699.05
5 0.44 50.0 125,000 143,765.83 140,000 408,765.83 419,999.99 11,234.16
6 0.45 50.2 125,000 147,116.59 140,000 412,116.59 421,679.99 9,563.4
7 0.46 51.4 125,000 150,476.95 140,000 415,476.95 431,759.99 16,283.04
8 0.46 50.22 125,000 151,206.63 140,000 416,206.63 421,847.99 5,641.36
9 0.46 54.16 125,000 151,456.26 140,000 416,456.26 454,943.99 38,487.74
10 0.48 54.21 125,000 156,650.41 140,000 421,650.41 455,363.99 33,713.58
11 0.48 49.78 125,000 156,813.63 140,000 421,813.63 418,151.99 -3,661.64
12 0.46 52.29 125,000 152,387.56 140,000 417,387.56 439,235.99 21,848.44
13 0.46 49.16 125,000 152,445.16 140,000 417,445.16 412,943.99 -4,501.17
14 0.46 54.44 125,000 151,869.1 140,000 416,869.1 457,295.99 40,426.89
15 0.46 52.85 125,000 151,091.42 140,000 416,091.42 443,939.99 27,848.58
16 0.46 49.37 125,000 151,513.86 140,000 416,513.86 414,707.99 -1,805.87
17 0.47 54.5 125,000 153,981.33 140,000 418,981.33 457,799.99 38,818.66
18 0.45 56.6 125,000 148,076.69 140,000 413,076.69 475,439.99 62,363.3
19 0.44 50.74 125,000 145,225.19 140,000 410,225.19 426,215.99 15,990.81
20 0.43 44.76 125,000 142,690.52 140,000 407,690.52 375,983.99 -31,706.52
21 0.42 44.1 125,000 138,658.08 140,000 403,658.08 370,439.99 -33,218.09
22 0.4 42.64 125,000 131,399.7 140,000 396,399.7 358,175.99 -38,223.71
23 0.4 36.88 125,000 130,804.44 140,000 395,804.44 309,791.99 -86,012.45
24 0.4 45.55 125,000 132,446.22 140,000 397,446.22 382,619.99 -14,826.22
25 0.4 43.77 125,000 130,948.45 140,000 395,948.45 367,667.99 -28,280.46
26 0.39 45.67 125,000 128,653.81 140,000 393,653.81 383,627.99 -10,025.82
27 0.39 45.66 125,000 128,692.21 140,000 393,692.21 383,543.99 -10,148.22
28 0.39 42.83 125,000 128,068.15 140,000 393,068.15 359,771.99 -33,296.15
29 0.38 40.99 125,000 124,314.14 140,000 389,314.14 344,315.99 -44,998.15
30 0.38 42.59 125,000 123,776.49 140,000 388,776.49 357,755.99 -31,020.49
In [105]:
print('Average daily profit = ',df_part1.Profit.mean().round(2),'GBP')
Average daily profit =  4126.22 GBP
In [106]:
print('Price volatility = ',df_part1.Profit.std().round(2),'GBP')
Price volatility =  36787.74 GBP
In [107]:
print('Profit for the whole period = ',df_part1.Profit.sum().round(2),'GBP')
Profit for the whole period =  123786.59 GBP

Part 2

In part 2, I will estimate the 95% expected shortfall for daily profit, based on the market prices for gas and electricity from the given dataset for November 2019. I will use bivariate normal distribution, with means and covariances taken from the dataset.

Preparing parameters

I will check for statistical parameters from a given dataset: means and covariances, to use it in creating bivariate normal distribution.

In [108]:
# Means for gas and electricity prices
means = (data.gas.mean(),data.electricity.mean())
print('Gas price mean:', means[0].round(2),'GBP')
print('Electricity price mean:', means[1].round(2),'GBP')
Gas price mean: 0.43 GBP
Electricity price mean: 49.0 GBP
In [109]:
# Covariance matrix of gas and electricity variables
z = np.array([data.gas, data.electricity])
covMatrix = np.cov(z, bias=True) # with switched bias function to normalize values
pd.DataFrame(covMatrix,columns=['Gas','Electricity'])
Out[109]:
Gas Electricity
0 0.000939 0.119763
1 0.119763 26.475282
In [110]:
print('Covariance (gas,electricity) = ',covMatrix[0][1])
Covariance (gas,electricity) =  0.11976307917849824

Simulation

I will generate 10'000 observations from bivariate normal distribution, with given parameters, to create the sample.

In [111]:
# Creating bivariate normal distribution with parameters from the dataset
mn = multivariate_normal(mean=means, cov=covMatrix)
In [112]:
# N - Number of observations
N = 10000

# Creating a sample of 10'000 observations from bivariate normal distribution with seed (random_state) for repeatable results
sample = mn.rvs(size=N, random_state=888)

pd.DataFrame(sample,columns=['Gas Price /therm [GBP]','Electricity Price /MWh [GBP]']).round(2)
Out[112]:
Gas Price /therm [GBP] Electricity Price /MWh [GBP]
0 0.43 49.91
1 0.41 44.75
2 0.45 52.36
3 0.42 47.88
4 0.45 48.49
... ... ...
9995 0.44 51.48
9996 0.46 52.32
9997 0.39 41.51
9998 0.41 48.98
9999 0.39 43.26

10000 rows × 2 columns

In [113]:
# Saving the sample for marker's review
pd.DataFrame(sample,columns=['Gas Price /therm [GBP]','Electricity Price /MWh [GBP]']).to_csv("sample_bivariate.csv")
In [114]:
# Visualization of the sample's distribution
sns.jointplot(x=sample[:,0],
              y=sample[:,1], 
              kind="kde", 
              space=0);

Calculating profit from created sample

In [115]:
# Creating dataframe from the sample
df = pd.DataFrame(sample,columns=['gas','electricity'])

# Gas consumption with contract price
df['Gas cost (contract)']=125000

# Gas consumption with market price
df['Gas cost (market)']= (578671.80-250000)*df.gas

# Operation daily cost
df['Operation cost']=140000

# Total Cost
df['Total Cost']=df['Gas cost (contract)']+df['Gas cost (market)']+df['Operation cost']

# Sales
df['Sales']=(8399999.85/1000)*df.electricity

# Profit
df['Profit']=df['Sales']-df['Total Cost']

# Loss - negative Profit for ES calculations
df['Loss']= -df['Profit']

df.head().round(2).style.format("{:,}")
Out[115]:
gas electricity Gas cost (contract) Gas cost (market) Operation cost Total Cost Sales Profit Loss
0 0.43 49.91 125,000 142,593.04 140,000 407,593.04 419,224.2 11,631.16 -11,631.16
1 0.41 44.75 125,000 136,369.93 140,000 401,369.93 375,875.14 -25,494.8 25,494.8
2 0.45 52.36 125,000 148,163.86 140,000 413,163.86 439,810.3 26,646.44 -26,646.44
3 0.42 47.88 125,000 136,969.57 140,000 401,969.57 402,196.23 226.66 -226.66
4 0.45 48.49 125,000 148,919.76 140,000 413,919.76 407,275.13 -6,644.63 6,644.63

Calculating VaR

For Value at Risk at alpha = 95% level, I will use Percent Point Function (PPF), inverse of Cumulative Distribution Function (CDF), from the sample's Loss values parameters - mean and standard deviation (std).

In [116]:
alpha = 0.95

VaR_2 = norm.ppf(alpha, np.mean(df.Loss), np.std(df.Loss))

print('Loss Value at Risk at 95% confidence level:', f"{round(VaR_2,2):,}",'GBP')
Loss Value at Risk at 95% confidence level: 56,582.15 GBP

Calculating Expected Shortfall

I will calculate 95% Expected Shortfall in two ways: discrete (empirical values from the sample Loss values) and continuous (from Probability Distribution Function, created with parameters from the sample Loss values). We will see that these two approaches give similar results, with difference due to normal distribution approximation.

In [117]:
# Discrete approach - average all Losses greater than 95% VaR

ES_empirical_2 = df[df.Loss>=VaR_2]['Loss'].mean() 

print('Expected Shortfall at 95% confidence level from sample:', f"{round(ES_empirical_2,2):,}",'GBP')
Expected Shortfall at 95% confidence level from sample: 71,908.99 GBP
In [118]:
# Continuous approach - calculating expected 

alpha = 0.05

ES_cont_2 =  np.mean(df.Loss) + np.std(df.Loss) * norm.pdf(norm.ppf(1-alpha))/alpha

print('Expected Shortfall at 95% confidence level from normal distribution function approximation:',f"{round(ES_cont_2,2):,}",'GBP')
Expected Shortfall at 95% confidence level from normal distribution function approximation: 71,842.21 GBP
In [119]:
# Visualization of sample Loss disitrbution
mean = np.mean(df.Loss)
std_dev = np.std(df.Loss)

df.Loss.hist(bins=1000, normed=True, histtype='stepfilled', alpha=0.5)

x = np.linspace(mean - 3*std_dev, mean + 3*std_dev, 100)

plt.plot(x,norm.pdf(x, mean, std_dev), "black" )
plt.title('Sample\'s Losses distribution')
plt.xlabel('Loss GBP')
plt.ylabel('Probaility')
plt.axvline(VaR_2,color='r')
plt.axvline(ES_empirical_2,color='g')
plt.text(VaR_2+3000,0.00001,'VaR 95 = '+str("{:,}".format(int(VaR_2))),rotation='vertical',color='r')
plt.text(ES_empirical_2+3000,0.00001,'ES 95 = '+str("{:,}".format(int(ES_empirical_2))),rotation='vertical',color='g')
plt.show()
C:\Users\Longin\Anaconda3\lib\site-packages\pandas\plotting\_matplotlib\hist.py:316: MatplotlibDeprecationWarning:


The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.

Conclusion:

The difference between empirical (ES_empirical) and continuous (ES_cont) Expected Shortfall figures is due to normal distribution approximations. However, we see that the Expected Shortfall Loss will be on average 71'843.97 GBP.

Part 3

In part 3 of the assignment, I use simulation approach, like in Part 2, to generate observations from bivariate normal distribution, with given mean and covariances parameters. Next, I modified 4% of all observations (400 days) to include generator failure (half production capacity). I choose fixed number of the sample failures, so variant B from the assignment guidelines. I have changed following parameters for failure days:

  • Gas consumption on a failure day = 578 671,80 / 2 = 289 335,90 therms / day
  • Electricity production on a failure day = 8400 MWh / 2 = 4200 MWh / day
In [120]:
# Generate 4% unique random numbers within a range 0 to N - each number correspond to a failure day in the sample
random.seed(888)
failure_days = random.sample(range(0, N), int(N/25))
In [121]:
# Copy of the sample from Part 2
df_with_failures = df.copy()
In [122]:
# Modifying sample to include 400 failure days (4% of total observations)

for i in failure_days:
    
    df_with_failures.at[i,'Gas cost (market)']=((578671.80*0.5)-250000)*df_with_failures.iloc[i]['gas']
    
    df_with_failures.at[i,'Total Cost']=df_with_failures.iloc[i]['Gas cost (contract)']+df_with_failures.iloc[i]['Gas cost (market)']+df_with_failures.iloc[i]['Operation cost']
    
    df_with_failures.at[i,'Sales']=(8399999.85/2/1000)*df_with_failures.iloc[i]['electricity']
    
    df_with_failures.at[i,'Profit']=df_with_failures.iloc[i]['Sales']-df_with_failures.iloc[i]['Total Cost']

df_with_failures['Loss']=-df_with_failures['Profit']

df_with_failures
Out[122]:
gas electricity Gas cost (contract) Gas cost (market) Operation cost Total Cost Sales Profit Loss
0 0.433846 49.907643 125000 142593.039910 140000 407593.039910 419224.197238 11631.157328 -11631.157328
1 0.414912 44.747041 125000 136369.932020 140000 401369.932020 375875.135947 -25494.796074 25494.796074
2 0.450796 52.358370 125000 148163.861688 140000 413163.861688 439810.299817 26646.438129 -26646.438129
3 0.416737 47.880505 125000 136969.573658 140000 401969.573658 402196.232513 226.658855 -226.658855
4 0.453096 48.485135 125000 148919.759315 140000 413919.759315 407275.130219 -6644.629096 6644.629096
... ... ... ... ... ... ... ... ... ...
9995 0.435036 51.481982 125000 142984.159943 140000 407984.159943 432448.640573 24464.480630 -24464.480630
9996 0.460554 52.315929 125000 151371.111273 140000 416371.111273 439453.797457 23082.686185 -23082.686185
9997 0.385268 41.506691 125000 126626.675743 140000 391626.675743 348656.198200 -42970.477543 42970.477543
9998 0.406562 48.979993 125000 133625.510196 140000 398625.510196 411431.930090 12806.419894 -12806.419894
9999 0.390512 43.256896 125000 128350.340466 140000 393350.340466 363357.922914 -29992.417552 29992.417552

10000 rows × 9 columns

Calculating VaR (as in part 2)

In [123]:
alpha = 0.95

VaR_3 = norm.ppf(alpha, np.mean(df_with_failures.Loss), np.std(df_with_failures.Loss))

print('Loss Value at Risk at 95% confidence level:', f"{round(VaR_3,2):,}",'GBP')
Loss Value at Risk at 95% confidence level: 64,388.32 GBP

Calculating Expected Shortfall (as in part 2)

In [124]:
# Discrete approach - average all Losses greater than 95% VaR

ES_empirical_3 = df_with_failures[df_with_failures.Loss>=VaR_3]['Loss'].mean() 

print('Expected Shortfall at 95% confidence level from sample:', f"{round(ES_empirical_3,2):,}",'GBP')
Expected Shortfall at 95% confidence level from sample: 81,690.71 GBP
In [125]:
# Continuous approach - calculating expected 

alpha = 0.05

ES_cont_3 =  np.mean(df_with_failures.Loss) + np.std(df_with_failures.Loss) * norm.pdf(norm.ppf(1-alpha))/alpha

print('Expected Shortfall at 95% confidence level from normal distribution function approximation:',f"{round(ES_cont_3,2):,}",'GBP')
Expected Shortfall at 95% confidence level from normal distribution function approximation: 80,816.11 GBP

The difference between empirical (ES_empirical) and continuous (ES_cont) Expected Shortfall figures is due to normal distribution approximations. However, we see that the Expected Shortfall Loss, for a sample with 4% of failure days, will by average 81'254.42 GBP.

In [126]:
mean = np.mean(df_with_failures.Profit)
std_dev = np.std(df_with_failures.Profit)

df_with_failures.Profit.hist(bins=1000, normed=True, histtype='stepfilled', alpha=0.5)

x = np.linspace(mean - 3*std_dev, mean + 3*std_dev, 100)

plt.plot(x, norm.pdf(x, mean, std_dev), "black" )
plt.title('Sample\'s Losses distribution (with failure days)')
plt.xlabel('Loss GBP')
plt.ylabel('Probaility')
plt.axvline(VaR_3,color='r')
plt.axvline(ES_empirical_3,color='g')
plt.text(VaR_3+3000,0.00001,'VaR = '+str("{:,}".format(int(VaR_3))),rotation='vertical',color='r')
plt.text(ES_empirical_3+3000,0.00001,'ES = '+str("{:,}".format(int(ES_empirical_3))),rotation='vertical',color='g')
plt.show()
C:\Users\Longin\Anaconda3\lib\site-packages\pandas\plotting\_matplotlib\hist.py:316: MatplotlibDeprecationWarning:


The 'normed' kwarg was deprecated in Matplotlib 2.1 and will be removed in 3.1. Use 'density' instead.

Conclusion:

We have obtained VaR 95% of 64,388.32 GBP and Expected Shortfall 95% of 81'254.42 GBP. The numbers makes sense in comparison to the sample without failure days - we expected increase in both VaR and ES in part 2, as failure days in the new sample generate more losses.

Part 4

In part 4 of the assignment, I will minimize Expected Shortfall, by hedging with Contracts for Differences for gas and electricity. To calculate optimal quantities for CFDs, I will search in feasible region (i.e. all allowed CFD quantities combinations) and pick quantities, which yield minimum Expected Shortfall value. I will use the grid search algorithm, presented below as pseudocode.

Algorithm

QG in range (0, max energy production 8`400 MWh )
QE in range (0, max consumption of gas (578`671.80 therms/day)

for x in in QG:
    for y in QE:
        calculate Expected Shortfall(x,y)

Find min(Expected Shortfall) and its corresponding (x,y) = (CFD_GAS_QTY,CFD_ELEC_QTY)

CFDs profits calculations

  • Profit from Gas CDF = (gas market price - CFD gas strike price) * gas quantity
  • Profit from Electricity CDF = (CFD electricity strike price - electricity market price) * electricity quantity
In [127]:
df_part4 = df_with_failures.copy()
df_part4['Gas Qty'] = 578671.80
df_part4['Elec Qty'] = 8400
df_part4.head().round(2).style.format("{:,}")
Out[127]:
gas electricity Gas cost (contract) Gas cost (market) Operation cost Total Cost Sales Profit Loss Gas Qty Elec Qty
0 0.43 49.91 125,000 142,593.04 140,000 407,593.04 419,224.2 11,631.16 -11,631.16 578,671.8 8,400
1 0.41 44.75 125,000 136,369.93 140,000 401,369.93 375,875.14 -25,494.8 25,494.8 578,671.8 8,400
2 0.45 52.36 125,000 148,163.86 140,000 413,163.86 439,810.3 26,646.44 -26,646.44 578,671.8 8,400
3 0.42 47.88 125,000 136,969.57 140,000 401,969.57 402,196.23 226.66 -226.66 578,671.8 8,400
4 0.45 48.49 125,000 148,919.76 140,000 413,919.76 407,275.13 -6,644.63 6,644.63 578,671.8 8,400
In [128]:
for i in failure_days:
    
    df_part4.at[i,'Elec Qty'] = 8400/2
    df_part4.at[i,'Gas Qty'] = 578671.80/2
    
df_part4['Total_Cost/MWh'] = df_part4['Total Cost']/df_part4['Elec Qty']
In [129]:
df_part4.head(5).round(2).style.format("{:,}")
Out[129]:
gas electricity Gas cost (contract) Gas cost (market) Operation cost Total Cost Sales Profit Loss Gas Qty Elec Qty Total_Cost/MWh
0 0.43 49.91 125,000 142,593.04 140,000 407,593.04 419,224.2 11,631.16 -11,631.16 578,671.8 8,400 48.52
1 0.41 44.75 125,000 136,369.93 140,000 401,369.93 375,875.14 -25,494.8 25,494.8 578,671.8 8,400 47.78
2 0.45 52.36 125,000 148,163.86 140,000 413,163.86 439,810.3 26,646.44 -26,646.44 578,671.8 8,400 49.19
3 0.42 47.88 125,000 136,969.57 140,000 401,969.57 402,196.23 226.66 -226.66 578,671.8 8,400 47.85
4 0.45 48.49 125,000 148,919.76 140,000 413,919.76 407,275.13 -6,644.63 6,644.63 578,671.8 8,400 49.28

Finding minimal value of ES

I used grid search approach, so looking for all possible combinations of quantities for Gas and Electricity contracts to find minimum Expected Shortage figure.

In [130]:
# Strike prices for CFDs
Q_strike_Gas = 0.43 # 0.43 GBP per gas therm
Q_strike_Elec = 49 # 49 GBP per MWh

# Feasible region in minimization problem

# Quantity range of CFD for gas - <0,587'671>
Q_G = range(0,587671,500)

# Quantity range of CFD for electricity - <0,8'400>
Q_E = range(0,8400,50)

# Confidence level for 95% Expected Shortfall 
alpha = 0.05 

# List for Expected Shortfall values for every CFDs combination
es = []

for qg in Q_G:
    for qe in Q_E:
        
        # Profit from CFD for Electricity with given quantity qe
        df_part4['Profit CFD Elec'] = (Q_strike_Elec - df_part4['electricity']) * qe
        
        # Profit from CFD for Gas with given quantity qg
        df_part4['Profit CFD Gas'] = (df_part4.gas - Q_strike_Gas) * qg
        
        df_part4['Total_Profit'] = df_part4.Profit + df_part4['Profit CFD Elec'] + df_part4['Profit CFD Gas']
        
        df_part4['Total_Loss'] = -df_part4['Total_Profit']
        
        es.append((norm.pdf(norm.ppf(alpha))* 1/alpha * np.std(df_part4.Total_Loss) - np.mean(df_part4.Total_Loss), qg, qe))
In [131]:
# Saving Expected Shortfall results to csv file
pd.DataFrame(es,columns=['ES','Qty_CFD_Gas','Qty_CFD_Electricity']).to_csv('min_es.csv')
In [132]:
df_es = pd.DataFrame(es,columns=['ES','Q_G','Q_E'])
In [133]:
min_ES = df_es[df_es.ES==df_es['ES'].min()].iloc[0][0].round(2)
CFD_Gas_Qty = df_es[df_es.ES==df_es['ES'].min()].iloc[0][1].round(2)
CFD_Electricity_Qty = df_es[df_es.ES==df_es['ES'].min()].iloc[0][2]
In [134]:
print('Minimum Expected Shortfall = ', f"{min_ES:,}",'GBP')
Minimum Expected Shortfall =  34,957.11 GBP
In [135]:
print('CFD Gas Qty for min Expected Shortfall = ', f"{CFD_Gas_Qty:,}",'gas therms')
CFD Gas Qty for min Expected Shortfall =  245,500.0 gas therms
In [136]:
print('CFD Electricity Qty for min Expected Shortfall = ', f"{CFD_Electricity_Qty:,}",'MWh')
CFD Electricity Qty for min Expected Shortfall =  7,900.0 MWh
In [137]:
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.tri as mtri

triang = mtri.Triangulation(df_es['Q_G'], df_es['Q_E'])

fig = plt.figure(figsize=(20,10))


ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_trisurf(triang, df_es['ES'], cmap='RdYlBu', linewidth=1)
ax.view_init(elev=20, azim=50)
ax.dist = 10


ax.set_xlabel('\n\n\nGas [therms]')
ax.set_ylabel('\n\n\nElectricity [MWh]')
ax.set_zlabel('\n\nExcpeted Shortfall [GBP]')

fig.colorbar(surf, shrink=0.5, aspect=10)

fig.suptitle('Optimisation for Expected Shortfall')
plt.show()

The below 3D graph generated in plotly, allows to rotate the figure and check where exactly the global minimum is. However, it can be only generated dynamically in Python or jupyter notebook, therefore I include a photo.

In [140]:
import plotly.graph_objects as go

min_ES_index = df_es[df_es.ES==df_es['ES'].min()].index[0]

fig = go.Figure(data=[go.Scatter3d(z=df_es['ES'], 
                                   x=df_es['Q_G'], 
                                   y=df_es['Q_E'],
                                   mode='markers', 
                                   marker=dict(size=1,  color=df_es["ES"])),
                     
                      # Min point
                      go.Scatter3d(z=[df_es['ES'][min_ES_index]], 
                                   x=[df_es['Q_G'][min_ES_index]], 
                                   y=[df_es['Q_E'][min_ES_index]], 
                                   mode='markers', 
                                   text='Minimum Expected Shortfall', 
                                   marker=dict(size=10))])

fig.update_layout(title='Optimization for Expected Shortfall', 
                  autosize=True,
                  width=800, height=800,
                  xaxis=dict(title_text="x"),
                  yaxis=dict(title_text="Electricity CFD Qty"))

fig.show()

Conclusion:

Minimizing 95% Expected Shortfall by checking all combinations of quantities for gas and electricity CFD contracts gave minimum Expected Shortfall of 34,957.11 GBP with accordingly 245,500 gas therms and 7,900 MWh electricity CFD's quantities. However, looking in the whole feasible region for the objective function was extremely computationally expensive, thus I had to constraint my search. I choose step of 500 for gas and 50 for electricity quantities, reducing computations and decreasing accuracy of the solution. Nevertheless, tested with different steps and ranges, the solution tend to descent to the global optimal solution, marked with the red dot on the 3d graph.

To summarize, in part 1 I calculated average daily profit over month for the generator of 4126.22 GBP. In part 2, I created a sample, drawing 10000 observations by random from bivariate normal distribution with parameters given in the dataset from November 2019. I obtained 56,582.15 GBP for VaR 95% and 71'843.97 GBP for Expected Shortfall 95%. In part 3, I introduced 400 observations to the sample from part 2, when only one generator unit is working, thus increasing overall losses risk. In fact, VaR 95% and Expected Shortfall 95% have increased in comparison to figures from sample not containing failure days, respectively up to 64,388.32 GBP and 81'254.42 GBP. In part 4, with the sample from part 3, I minimized Expected Shortfall 95% by searching for most optimal quantities of CFD contracts for gas and electricity. I obtained 34,957.11 GBP of Expected Shortfall 95%, largely reducing the losses risk, in comparison to sample without futures contracts, by around 57% (46,300 GBP). Therefore, I conclude that using CFD contract is be beneficial to the generator operator, reducing the losses risk heavily.